;;########################################################################
;; OLS, Robust & Monotonic Regression ViSta model object
;; Contains code for robust regression object
;; and for simple regression visualization
;; Copyright (c) 1995-19966 by Carla M. Bann
;; Copyright (c) 1997-2000 by Carla M. Bann and Forrest W. Young
;;########################################################################

;IN ORIGINAL ROBUST.LSP ... LOADED BEFORE ANYTHING EXCEPT MOB

(defproto ROBUST-PROTO
  '(weight-list
    c)
  ()
  morals-proto
  "Robust Regression Model")


(defmeth ROBUST-PROTO :weight-list (&optional (obj-id nil set))
  (when set (setf (slot-value 'weight-list) obj-id))
  (slot-value 'weight-list))

(defmeth ROBUST-PROTO :c (&optional (value nil set))
  (when set (setf (slot-value 'c) value))
  (slot-value 'c))


(defun robust-regression-model (x y &key 
                                   (intercept T)
                                   (print T)
                                   (plot T)
                                   weights
                                   (included (repeat t (length y)))
                                   predictor-names
                                   response-name
                                   case-labels                                  
                                   (iterations 0)
                                   vista-object
                                   (c 4.685)
                                   )
    (let* ((m (send ROBUST-PROTO :new x y)))
    (send m :x x)
    (send m :y y)
    (send m :intercept intercept)
    (send m :weights weights)
    (send m :included included)
    (send m :predictor-names predictor-names) 
    (send m :response-name response-name)
    (send m :case-labels case-labels)
    (send m :iterations iterations)
    (send m :vista-object vista-object)
    (send m :c c)
    (send m :compute)
    (send m :count 0)
    (setf *Robust-Model* m)
    m))

; :ISNEW method

(defmeth ROBUST-PROTO :isnew (x y)
  (call-method regression-model-proto :isnew)
  self)

 

;FROM TOP PORTION OF REGVITER.LSP 


(defmeth robust-proto :biweight (x)
  (^ (- 1 (^ (pmin (abs x) 1) 2)) 2))


(defmeth robust-proto :robust-weights (&optional (c 4.685))
  (let* ((vista-obj (send self :vista-object))
         (rr (send self :raw-residuals))
         (s  (/ (median (abs rr)) .6745)))
  (send self :biweight (/ rr (* c s)))))


(defmeth robust-proto :max-relative-error (x y)
  (max (/ (abs (- x y))
          (pmax (sqrt machine-epsilon)
                (abs y)))))

(defmeth robust-proto :robust-loop (&optional (epsilon .001)
                                              (limit 20)
                                              (c 4.685))
  (let* ((count 0)
        (last-beta nil)
        (beta (send self :coef-estimates))
        (rel-err 0)
        (weight-list nil))
  (if (equalp (send self :count) 0)
        (setf weight-list (list (send self :weights)))
        (setf weight-list (send self :weight-list)))
  (if (equalp (send self :count) 0)
      (send self :rsq-list (list (send self :r-squared)))) 
    (loop
     (send self :weights (send self :robust-weights c))
     (setf weight-list (remove 'nil (append weight-list (list (send self :weights)))))
     (send self :weight-list weight-list)
     (send self :count (+ (send self :count) 1))
     (setf count (+ count 1))
     (setf last-beta beta)
     (setf beta (send self :coef-estimates))
     (setf rel-err 
           (send self :max-relative-error beta last-beta))
     (send self :RSq-list (append (send self :RSq-list) 
                                    (list (send self :R-squared))))
     (when (or (< rel-err epsilon)
             (= limit count))
           (send self :compute)
           (return
            (list beta rel-err count))))))


(defmeth robust-proto :robust-bootstrap (&optional
                                         (nb 100)
                                         (epsilon .01))
  (let* ((n (send self :num-cases))
         (k (- (send self :num-coefs) 1))
         (x (send self :x))
         (y (send self :y))
         (result nil)
         (i-n (iseq n))
         (i-k (iseq k)))
    (dotimes (i nb)
       (let ((s (sample i-n n t)))
         (send self :x (select x s i-k))
         (send self :y (select y s))
         (send self :weights nil)
         (push (first (send self :robust-loop epsilon))
               result)))
    (send self :x x)
    (send self :y y)
    (send self :weights nil)
    (transpose result)))

